Grey radiation

This notebook describes the grey radiation model and goes through some examples of things that can be done with it.

The basis of this model is to produce the vertically one dimensional structure of the atmosphere given a concentration distribution of greenhouse gases.

I used chapters 10 and 14 of Atmospheric Circulation Dynamics and General Circulation Models and chapter 4 of Principles of Planetary Climate to help produce the model.

Model description

Conservation of energy

The total energy of the atmosphere is given by: \begin{equation} \rho e^{tot} = \frac{1}{2}\rho \pmb{v}^2 + \rho \Phi + \rho u \label{eq: 1}\tag{1} \end{equation} This is the sum of kinetic energy, gravitational potential energy ($\Phi = gz$) and internal energy. There is no source of total energy, meaning it is conserved: \begin{equation} \frac{\partial}{\partial t}(\rho e^{etot}) + \nabla \cdot \pmb{F}^{etot} = 0 \label{eq: 2}\tag{2} \end{equation} Where $\textbf{F}^{etot}$ is the flux density vector of total energy determined by the the balance equations for each component of $e^{tot}$.

\begin{equation} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \pmb{v}) = 0 \label{eq: 3}\tag{3} \end{equation}\begin{equation} \frac{d\pmb{v}}{dt} + 2\pmb{\Omega} \times \pmb{v} = -\frac{1}{\rho}\nabla p - \nabla \Phi \label{eq: 4}\tag{4} \end{equation}

Using $(2)$ along with conservation of mass $(3)$ and conservation of momentum $(4)$ as well as the fact $\frac{\partial \Phi}{\partial t}=0$, we can obtain an equation for the evolution of kinetic and potential energy in flux form (note that $\frac{d}{dt} = \frac{\partial}{\partial t} + \pmb{v} \cdot \nabla$):

\begin{equation} \frac{\partial}{\partial t}(\frac{1}{2}\rho \pmb{v}^2) + \nabla \cdot \left[ \left(\frac{1}{2}\rho \pmb{v}^2 + p\right)\pmb{v} \right] = p \nabla \cdot \pmb{v} - \rho \cdot \nabla \Phi \label{eq: 5}\tag{5} \end{equation}\begin{equation} \frac{\partial}{\partial t}(\rho \Phi) + \nabla \cdot \left(\rho \Phi \pmb{v}\right) = \rho \pmb{v}\cdot \nabla \Phi \label{eq: 6}\tag{6} \end{equation}

We can determine the internal energy equation by noting that the right hand sides (RHS) of $(5)$ and $(6)$ are source terms of energy so the source term for the internal energy equation must be $-p \nabla \cdot \pmb{v}$ so the source term of total energy (sum of RHS of $(5)$, $(6)$ and $(7)$) equals $0$ as to satisfy equation $(2)$. Thus the internal energy equation is:

\begin{equation} \frac{\partial}{\partial t}(\rho u) + \nabla \cdot \left(\rho u \pmb{v} + \pmb{F}\right) = -p \nabla \cdot \pmb{v} \label{eq: 7}\tag{7} \end{equation}

$\pmb{F}$ is the flux density of internal energy other than the advection term i.e. it is the heat flux. In the grey radiation model, we assume the only contribution is from radiation (as opposed to sensible heat - molecular thermal diffusion).

We want to use all these equations to get an expression for how temperature, $T$, responds to changes in $\pmb{F}$. We can combine $(7)$ and $(3)$ to give:

\begin{equation} \rho \left(\frac{du}{dt} + p\frac{dv_s}{dt}\right) = -\nabla \cdot \pmb{F} \label{eq: 8}\tag{8} \end{equation}

Where $v_s = 1/ \rho$ is specific volume. If we then use the definition of enthalpy, $dh = du + d(pv_s)$, we get an expression for the evolution of enthalpy:

\begin{equation} \rho \frac{dh}{dt} = \frac{dp}{dt} - \nabla \cdot \pmb{F} \label{eq: 9}\tag{9} \end{equation}

We can then use $h = c_pT$ where $c_p$ is the heat capacity of air at constant pressure per unit mass (units are J/K/kg). We are also only interested in a 1D vertical ($\nabla \rightarrow \frac{\partial}{\partial z}$), static ($\pmb{v} \rightarrow 0$) atmosphere in this section. Under these conditions, $\frac{d}{dt} \rightarrow \frac{\partial}{\partial t}$, $(3)$ gives $\rho = constant$, $(4)$ gives $\frac{\partial p}{\partial t} = 0$ and $\frac{\partial}{\partial z} = -g\rho \frac{\partial}{\partial p}$. This gives our final temperature equation:

\begin{equation} \frac{\partial T}{\partial t} = \frac{g}{c_p}\frac{\partial F}{\partial p} \label{eq: 10}\tag{10} \end{equation}

In the model, we discretize this using a finite volume method (pressure and flux are at the boundaries of a cell whereas temperature is at the cell centre):

\begin{equation} \frac{T^{n+1}_{i+1/2} - T^{n}_{i+1/2}}{\Delta t} = \frac{g}{c_p}\frac{F^{n}_{i+1} - F^{n}_i}{p_{i+1} - p_i} \label{eq: 11}\tag{11} \end{equation}

Where $i$ indicates grid position and $n$ indicates time.

Flux equation

Now to complete the equation, we want to find an expression for $F$ as a function of $T$ and $p$.

Lets consider a slab of atmosphere of thickness $\Delta z$. Let there be an upwelling flux, $F^{\uparrow}_{\nu}(z)$ at the bottom of the slab and an upwelling flux $F^{\uparrow}_{\nu}(z+\Delta z)$ at the top. There is a change in flux because the layer of atmosphere absorbs and emits radiation. According to Kirchoff's law, emissivity = absorptivity = $e_{\nu}$ and we have:

\begin{equation} F^{\uparrow}_{\nu}(z+\Delta z) = (1-e_{\nu})F^{\uparrow}_{\nu}(z) + e_{\nu}\pi B(\nu, T) \label{eq: 12}\tag{12} \end{equation}

This means that for each frequency $\nu$, the radiation flux that reaches the top of this slab of atmosphere is equal to the amount that was incident upon the bottom of the slab minus the amount absorbed by the slab plus the amount emitted by the slab. We assume local thermodynamic equilibrium so the emitted radiation follows the blackbody distribution, $B(\nu, T)$.

The emissivity is given by $e_{\nu} = \rho (z)k(\nu, z)\Delta z$. $\rho$ is the density of the fluid ($kgm^{-3}$). $k$ is an absorption coefficient and has units of area per unit mass ($m^2kg^{-1}$). It can be thought of as the area taken out of the incident beam by the absorbers contained in a unit mass of atmosphere. For a fluid with n different absorbing substances, $k = \sum_{i=0}^n k_iq_i$ where $q_i(z)$ is the specific humidity ($q_i = \rho_i / \rho$).

If we sub this into $(12)$, and take the limit $\Delta z \rightarrow 0$:

\begin{equation} \frac{dF^{\uparrow}_{\nu}}{dz} = -\rho kF^{\uparrow}_{\nu} + \rho k\pi B(\nu, T) \label{eq: 13}\tag{13} \end{equation}

If we now define optical depth, $\tau_{\nu}(z)$ such that $d\tau_{\nu} = \rho kdz$, we get:

\begin{equation} \frac{dF^{\uparrow}_{\nu}}{d\tau_{\nu}} = -F^{\uparrow}_{\nu} + \pi B(\nu, T) \label{eq: 14}\tag{14} \end{equation}

Optical depth increases as you move down through the atmosphere, coming from the surface ($\tau_{\nu} = 0$ at the surface and $\tau = \tau^{\infty}$ at the top of the atmosphere), and can be thought of as the effective distance travelled by light at a particular frequency given the interaction with the gas. The more the gas interacts with the light, the greater the optical depth. Going through the same analysis but considering the downwelling radiation incident upon the top of the slab, we get:

\begin{equation} F^{\downarrow}_{\nu}(z-\Delta z) = (1-e_{\nu})F^{\downarrow}_{\nu}(z) + e_{\nu}\pi B(\nu, T) \label{eq: 15}\tag{15} \end{equation}

which becomes:

\begin{equation} \frac{dF^{\downarrow}_{\nu}}{d\tau_{\nu}} = F^{\downarrow}_{\nu} - \pi B(\nu, T) \label{eq: 16}\tag{16} \end{equation}

However, because we have boundary conditions on $F^{\uparrow}$ and $F^{\downarrow}$ at the top of the atmosphere, it is convenient to redefine $\tau$ so that it is zero at the top of the atmosphere. We do this through: $\tau = \tau^{\infty} - \tau$. Subbing this into $(14)$ and $(16)$ gives us:

\begin{equation} \frac{dF^{\uparrow}_{\nu}}{d\tau_{\nu}} = F^{\uparrow}_{\nu} - \pi B(\nu, T) \label{eq: 17}\tag{17} \end{equation}\begin{equation} \frac{dF^{\downarrow}_{\nu}}{d\tau_{\nu}} = -F^{\downarrow}_{\nu} + \pi B(\nu, T) \label{eq: 18}\tag{18} \end{equation}

Grey gas approximation

The flux we need to sub into $(11)$ is $F = \int F^{\uparrow}_{\nu} d\nu - \int F^{\downarrow}_{\nu} d\nu$. So we need to know the frequency dependence of $\tau_{\nu}$. In the grey gas approximation, we assume that $k_{\nu}$ and thus $\tau_{\nu}$ doesn't depend on frequency. Integrating $(17)$ and $(18)$ over all frequencies then gives:

\begin{equation} \frac{dF^{\uparrow}}{d\tau} = F^{\uparrow} - \sigma_B T^4 \label{eq: 19}\tag{19} \end{equation}\begin{equation} \frac{dF^{\downarrow}}{d\tau} = -F^{\downarrow} + \sigma_B T^4 \label{eq: 20}\tag{20} \end{equation}

where $\sigma_B$ is the Boltzmann constant. The solution to $(19)$ is given by:

\begin{equation} F^{\uparrow}(\tau) = F^{\uparrow}(\tau = 0)e^{\tau} + \int_0^{\tau} \sigma_B T(\tau')^4e^{\tau - \tau'} d\tau' \label{eq: 21}\tag{21} \end{equation}

In the model, we find the flux at altitude level $i-1$, corresponding to $p + \Delta p$ and $\tau +\Delta \tau$ from the flux at paltitude level $i$, corresponding to $p$ and $\tau$. The equation for this is given by:

\begin{equation} F^{\uparrow}_{i} = F^{\uparrow}_{i-1}e^{-\Delta \tau} - \int_{\tau}^{\tau + \Delta \tau} \sigma_B T(\tau')^4e^{\tau - \tau'} d\tau' \label{eq: 22}\tag{22} \end{equation}

We simplify this, by assuming that $T$ is constant between $\tau$ and $\tau + \Delta \tau$, so we take it out the integral and give it the value at $i-\frac{1}{2}$. So we finally get:

\begin{equation} F^{\uparrow}_{i-1} = F^{\uparrow}_{i}e^{\Delta \tau} + \sigma_B T_{i-1/2}^4(1 - e^{\Delta \tau}) \label{eq: 23}\tag{23} \end{equation}

And similarly, for downwelling:

\begin{equation} F^{\downarrow}_{i-1} = F^{\downarrow}_{i}e^{-\Delta \tau} + \sigma_B T_{i-1/2}^4(1 - e^{-\Delta \tau}) \label{eq: 24}\tag{24} \end{equation}

We can then compute $F_i = F^{\uparrow}_{i} - F^{\downarrow}_{i}$ and sub into $(11)$.

The boundary conditions are $F^{\uparrow}(\tau = 0) = (1-A)F^{\odot}/4$ and $F^{\downarrow}(\tau = 0) = 0$. $F^{\odot} = 1367 Wm^{-2}$ is the solar constant and $A$ is the albedo, indicating how much of the solar radiation is reflected, a typical value is 0.3. The first boundary conditions basically states that the radiation emitted by the planet into space must balance the radiation incoming into the planet from the star. The factor of 4 is because the sunlight only impinges on an area $\pi R_{Earth}^2$ whereas the planet radiates from its entire area $4\pi R_{Earth}^2$. The second boundary condition states that the star does not provide any radiation in the frequency range emitted by the planet (planet is cooler than star so frequency spectrum is at much lower frequency). The fact that both boundary conditions are at the top of the atmosphere is the reason why equations $(23)$ and $(24)$ are written such that the flux at level i-1 is calculated from the flux at level i (not sure this makes physical sense for the upward flux though!!).

Short wavelength interacting atmosphere

Equations $(19)$ to $(24)$ are derived assuming the atmosphere only interacts with the long wavelength radiation emitted by the planet and by the atmosphere itself and not with the short wave radiation coming from the star.

To a first approximation, we can include the interaction with short wave radiation by instead of integrating $(17)$ and $(18)$ over all frequencies, but over two separate frequency regimes.

image.png

We get the short wave ($sw$) fluxes by integrating over the solar spectrum, so from about $\nu(\lambda = 0.1 \mu m)$ to $\nu(\lambda = 3 \mu m)$ as indicated above. We get the long wave ($lw$) fluxes by integrating over the eart spectrum, so from about $\nu(\lambda = 3 \mu m)$ to $\nu(\lambda = 100 \mu m)$. To do this we have to integrate over $B(\nu, T)$ in the two regimes. But the temperature in $B(\nu, T)$ is the temperature of a thin layer of the atmosphere and we expect this to be approximately the surface temperature of the planet and much smaller than the temperature of the star. Hence we can say:

\begin{equation} \int_{\nu(\lambda = 3 \mu m)}^{\infty} B(\nu, T) d\nu \approx 0 \label{eq: 25}\tag{25} \end{equation}\begin{equation} \int_{0}^{\nu(\lambda = 3 \mu m)} B(\nu, T) d\nu \approx \int_{0}^{\infty} B(\nu, T) d\nu = \sigma_B T^4 \label{eq: 26}\tag{26} \end{equation}

And our four equations are:

\begin{equation} \frac{dF^{\uparrow}_{lw}}{d\tau_{lw}} = F^{\uparrow}_{lw} - \sigma_B T^4 \label{eq: 27}\tag{27} \end{equation}\begin{equation} \frac{dF^{\downarrow}_{lw}}{d\tau_{lw}} = -F^{\downarrow}_{lw} + \sigma_B T^4 \label{eq: 28}\tag{28} \end{equation}\begin{equation} \frac{dF^{\uparrow}_{sw}}{d\tau_{sw}} = F^{\uparrow}_{sw} \label{eq: 29}\tag{29} \end{equation}\begin{equation} \frac{dF^{\downarrow}_{sw}}{d\tau_{sw}} = -F^{\downarrow}_{sw} \label{eq: 30}\tag{30} \end{equation}

The boundary conditions for the short wave equations are $F^{\uparrow}_{sw}(\tau = \tau_s) = AF^{\downarrow}_{sw}(\tau = \tau_s) = AF^{\odot}e^{-\tau_s}/4$ and $F^{\downarrow}_{sw}(\tau = 0) = F^{\odot}/4$. This is to satisfy the amount of solar light reflected by the planet and the flux of stellar light that reaches the planet respectively. Using these, we get:

\begin{equation} F^{\uparrow}_{sw} = \frac{1}{4}AF^{\odot}e^{\tau_{sw}-2\tau_{s}} \label{eq: 31}\tag{31} \end{equation}\begin{equation} F^{\downarrow}_{sw} = \frac{1}{4}F^{\odot}e^{-\tau_{sw}} \label{eq: 32}\tag{32} \end{equation}

The $F$ value to sub into $(11)$ in this case is $F = F^{\uparrow}_{lw} + F^{\uparrow}_{sw} - F^{\downarrow}_{lw} - F^{\downarrow}_{sw}$.

The key here is that $(27)-(28)$ and $(29)-(30)$ depend on different optical depths. So if an atmosphere had no gas that interacted with stellar light, $\tau_{sw}$ would be zero everywhere and $F_{sw}$ would be constant everywhere and just provide a boundary condition for $F_{lw}$ as in equations $(19)-(24)$ i.e. $\frac{\partial F_{sw}}{\partial p} = 0$ so $F_{sw}$ does not appear in equation $(10)$. But if an atmosphere did have a gas that interacted with stellar light, $\tau_{sw}$ would be positive at the surface and thus from $(32)$, it is clear that less radiation will reach the surface causing a cooling effect there.

Optical depth profiles

As in equation $(11)$, we have the flux at pressure level $i$, $F_i$, and all the flux equations are in terms of optical depth, $\tau$, we need a means of relating $\tau$ to pressure. The model gives four options in the file grey_optical_depth.py:

scale_height

This is probably the most physical and is the form used in Atmospheric Circulation Dynamics and General Circulation Models. The form of the function is:

\begin{equation} \tau = \tau_s \left(\frac{p}{p_s}\right)^{\alpha + 1} \label{eq: 33}\tag{33} \end{equation}

Where $\tau_s$ is the optical depth at the surface and $p_s$ is the pressure at the surface which is about 101320 Pa for Earth. The parameter $\alpha$, determines how sharply $\tau$ rises near the surface, a larger value means a steeper rise. $\alpha$ can be explained in more physical terms by considering the corresponding mass concentration of absorbing constituent profile:

\begin{equation} q = q_s \left(\frac{p}{p_s}\right)^{\alpha} = q_s e^{-z/h} \label{eq: 34}\tag{34} \end{equation}

Where $q_s$ is the value at the surface and the $q_s e^{-z/h}$ is valid in an isothermal atmosphere, in which case $h = H / \alpha$ where $h$ is the scale height of the mass concentration and $H$ is the pressure scale height.

The function is called using OpticalDepthFunctions.scale_height(p, p_width, tau_surface)[1], where the first argument is the pressure levels where you would like $\tau$ to be computed. p_width is the difference between the pressure level at the surface and the pressure level where mass concentration, $q$ falls to $q_s/e$. $\alpha$ is then computed from p_width. Note that all the optical depth functions have multiple outputs, hence we only take the second output here.

Some examples with changing p_width and tau_surface = 4 is shown on the left in the cell below while examples with changing tau_surface and p_width = 10000 are shown on the right.

exponential

The form of this is:

\begin{equation} \tau = \tau_s \frac{e^{\alpha p} - 1}{e^{\alpha p_s} - 1} \label{eq: 35}\tag{35} \end{equation}

It has the same arguments as [scale_height], so it is called through OpticalDepthFunctions.exponential(p, p_width, tau_surface)[1]. It gives a very similar profile, as shown in the cell below for p_width = 10000 and tau_surface = 4. The reason for including this function is explained in the equilibrium solutions section.

peak_in_atmosphere

The form of this is:

\begin{equation} \tau= \begin{cases} \tau_s \frac{e^{\alpha(p-p_{max})} - e^{-\alpha p_{max}}}{2 - e^{-\alpha p_{max}} - e^{\alpha (p_{max}-p_s)}} & p\leq p_{max}\\ \tau_s \frac{2 - e^{-\alpha p_{max}} - e^{\alpha (p_{max}-p)}}{2 - e^{-\alpha p_{max}} - e^{\alpha (p_{max}-p_s)}} & p > p_{max} \end{cases} \label{eq: 36}\tag{36} \end{equation}

This looks complicated, but it is derived from a simple mass concentration distribution that is peaked somewhere other than the surface:

\begin{equation} q = q(p_{max})e^{-\alpha |p-p_{max}|} \label{eq: 37}\tag{37} \end{equation}

The function is called through OpticalDepthFunctions.peak_in_atmosphere(p, p_width, p_max, tau_surface)[1]. p_max is the pressure level where the mass concentration distribution is peaked and now p_width is the difference between p_max and the pressure level where mass concentration, $q$ falls to $q(p_{max})/e$.

The cell below shows both the $\tau$ and $q$ distributions for p_max = 50000, tau_surface = 4 and varying p_width. Specific humidity, $q$ (see just above equation $(13)$) was derived assuming that the absoprtion coefficient, $k=1 m^2kg^{-1}$.

scale_height_and_peak_in_atmosphere

This is just the addition of scale_height with peak_in_atmosphere. It is called through OpticalDepthFunctions.scale_height_and_peak_in_atmosphere(p, p_width1, tau_surface1, p_width2, p_max2, tau_surface2)[1], where the arguments ending 1 correspond to scale_height and those ending 2 correspond to peak_in_atmosphere.

The below cell shows the $\tau$ and $q$ distributions for p_width1 = 1000, tau_surface1 = 4, p_width2 = 5000, p_max2 = 50000, tau_surface2 = 1.

Equilibrium analytic solutions

Equilibrium is reached when $\frac{\partial T}{\partial t} = 0$ so from $(10)$, we require $\frac{\partial F}{\partial p} = 0$ and hence $F$ is a constant and in the absence of heat supplied by the Earth through e.g. radioactive decay, we take the constant to be zero. So we have the four equations $(27) - (30)$ aswell as $F = F^{\uparrow}_{lw} + F^{\uparrow}_{sw} - F^{\downarrow}_{lw} - F^{\downarrow}_{sw} = 0$. That gives us 5 equations in 5 unknowns which we can solve.

Non short wave interacting

In the absence of any short wave interacting gas, $\tau_{sw} = 0$ and it becomes a lot simpler. The solution is below, where $F_{sw, eqb} = F^{\downarrow}_{sw, eqb} - F^{\uparrow}_{sw, eqb} = (1-A)F^{\odot}/4$: \begin{equation} F^{\uparrow}_{lw, eqb} = \frac{1}{2}F_{sw, eqb}(2 + \tau_{lw})\\ F^{\downarrow}_{lw, eqb} = \frac{1}{2}F_{sw, eqb}\tau_{lw}\\ T_{eqb} = \left[\frac{F_{sw, eqb}}{2\sigma_B}(1+\tau_{lw})\right]^{\frac{1}{4}} \label{eq: 38}\tag{38} \end{equation}

The below cell shows how you can view the analytic solution. First you need to create an instance of the GreyGas class by providing the optical depth function aswell as the required arguments. You also need to provide the number of grid points in pressure space, nz. If this is set to auto, the number of grid points will be chosen so as to get good coverage in both pressure and optical depth space through the function get_p_and_tau_interface_grids.

The equilibrium solution is then obtained through the function equilibrium_sol and is plotted through plot_eqb.

In the plot on the right, the negative fluxes correspond to $F^{\downarrow}$ while the positive fluxes correspond to $F^{\uparrow}$.

With short wave interacting

The solution to the 5 equations is a lot more complicated if $\tau_{sw} \neq 0$. First, you need to express all the variables in terms of one of the parameters $p$, $\tau_{lw}$ or $\tau_{sw}$. It is easiest if we use $\tau_{sw}$, and we end up having to compute the following integral, where $A$ is the albedo:

\begin{equation} \int_{0}^{\tau_{sw}} \frac{d\tau_{lw}}{d\tau_{sw}'} (e^{-\tau_{sw}'} - Ae^{\tau_{sw}'}) d\tau_{sw}' \label{eq: 39}\tag{39} \end{equation}

To get a solution to this, we require $\frac{d\tau_{lw}}{d\tau_{sw}'}$ is of the form $\sum_n (\tau_{sw}')^n$ where n is always an integer. The only way we can do this using our OpticalDepthFunctions is if both $\tau_{lw}$ and $\tau_{sw}$ are exponential and the ratio $\frac{\alpha_{lw}}{\alpha_{sw}}$ is a small integer (< 10, theoretically a solution is possible for larger integers but it takes a long time to compute).

If the optical depth functions supplied do not satisfy this requirement, the returned solution will be that with $\tau_{sw} = 0$ and the value of correct_solution will be set to False to indicated the wrong equilibrium has been calculated.

The below cell shows the equilibrium solution obtained by adding a short wave component to the lw_gas in the above cell. To ensure that the $\alpha$ ratio is an integer, we first need to find $\alpha_{lw}$ using get_exponential_alpha, then divide this by an integer before finding the corresponding $p_{width}$ using get_exponential_p_width.

In these plots, the dotted lines show the solution, where $\tau_{sw} = 0$.

It can be seen from this plot that the effect of adding the short wave component is to decrease the temperature at the surface and increase the temperature in the higher atmosphere. The surface temperature decreases because the short wave gas absorbs sunlight, reducing the amount that reaches the surface, as shown by the downward $F_{sw}$ in the plot on the right. This explains why the temperature of the sea decreases as you go deeper. The increase in temperature higher up is because the upward flux of short wave radiation has increased.

Physically, you can think that the short wave gas is heated by sunlight so if the short wave gas is dominant, the temperature will increase as you move to the top of the atmosphere. Long wave gas is heated by terrestrial light so if long wave gas is dominant, the temperature will increase as you move towards the surface.

Evolution to equilibrium

We can use equation $(11)$ to investigate how temperature evolves with time if the initial state of the atmosphere was not in equilibrium.

Long wave only

The function to call to do this is evolve_to_equilibrium and it returns a dictionary, dict, where dict['t'] gives a list of all the times of each step of the simulation and dict['T'][i] gives the temperature profile at time dict['t'][i].

An example for lw_gas initialised in the previous section is shown below. The temperature is adjusted at each time step using the function take_time_step(t, T_initial). t is just the current time and T_initial is the starting temperature profile, and is only necessary if t=0. If T_initial is not provided, the initial profile will be an isothermal atmosphere such that $F^{\uparrow}_{lw} + F^{\uparrow}_{sw} = F^{\downarrow}_{sw}$ at the top of the atmosphere.

This function will find the fluxes corresponding to the current temperature profile and then automatically find a time step such that exactly one pressure level has a temperature change of 1 Kelvin and all other levels have a smaller temperature change. This time step is found through the function update_time_step and the max temperature change will be lowered from 1 Kelvin if convergence isn't happening.

take_time_step returns the time updated by the timestep aswell as delta_net_flux which indicates by how much the net flux, $F$, is changing. The idea being that when this value is small equilibrium has been reached. In fact, the functioncheck_equilibrium(delta_net_flux, flux_thresh) explicitly checks for equilibrium by returning True if delta_net_flux < flux_thresh.

The evolution of the temperature is plotted by getting an array of all the times and temperatures then running Animate(grey_world, T_array, t_array, T_eqb, correct_solution, log_axis=True, nPlotFrames=100) using the class in the animation.py file. This then finds nPlotFrames (100 by default) suitable times to plot from all those in T_array. The log_axis parameter determines whether to plot pressure on a log scale, it is set to True by default. T_eqb does not have to be supplied but helps check if the analytical solution is being approached.

Short wave analytic

We can do the same thing for the sw_gas initialised just below equation $(36)$, as given by the next cell, for which we also have a correct analytical solution.

Stratosphere

If we don't require an analytic solution, we can obtain more interesting atmospheric structures. In the cell below, we include a gas that absorbs short wave radiation at a pressure level of 2,000 Pa. This then causes a temperature inversion, analagous to the role of ozone in our stratosphere.

Note that in the animation, the red curve doesn't converge to the orange curve because the orange solution is that with $\tau_{sw} = 0$ as indicated by the warning.

Mesosphere

The next example shows that if we include a peak in the mass concentration of the long wave gas above the peak of the short wave gas, we can construct an atmosphere which has a mesosphere where there is another temperature inversion.

In this example, we also track the fluxes at each step of the simulation. We do this by initialising the dictionary to pass to evolve_to_equilibrium. We must provide the first entry for each key of the dictionary. The data corresponding to the flux key is also a dictionary containing the up/down long-wave/short-wave permutations.

Thermosphere

Finally, if we include a short wave gas which has pretty much constant non-zero concentration except for a peak at a level above the mesosphere peak in long wave gas concentration, we can get a thermosphere, where the temperature increases as you move out to space.

Convective adjustment

We can also include convective adjustment which ensures that $dT/dz \geq -\gamma$ where $\gamma$ is the critical lapse rate, which by default is the dry adiabatic lapse rate, $\gamma = g/c_p$.

This is done by setting convective_adjust=True in evolve_to_equilibrium, take_time_step or equilibrium_sol.

The below shows a comparison between the final temperature profile for the thermosphere example above when convection is and is not included. You can see that in the locations where the temperature is increasing with increasing with increasing pressure (around $10^3 Pa$ and above $2 \times 10^4 Pa$), the rate of increase in the convective case is reduced.

The convective adjustment algorithm is explained in Convective Adjustment.ipynb.

Evolution with changing optical depth profile

We can also see how the Temperature profile changes in response to changes in the optical depth profile.

Example 1 - aerosols in stratosphere

In this example, we start off with just a long wave gas. We then increase the surface optical thickness with time until it reaches 6. At this stage we then let the atmosphere reach equilibrium. Once equilibrium is reached, we introduce a short wave gas peaked at 20000 Pa. We increase the short wave surface optical depth with time until it reaches 1.2. Then we let the atmosphere reach its new equilibrium. Finally, we remove all the short wave gas so the optical depth goes to 0 instantly and see the final equilibrium reached.

The idea behind this simulation is the geoengineering solution to global warming of inject aerosols into the stratosphere. So the initial stage indicates global warming where more greenhouse gases increase the surface temperature. The next step is the injection of aerosols which interact with short wave radiation and we can see from the animation below that the surface temperature starts decreasing when the short wave optical depth is non-zero. The final removal of short wave optical depth is included because aerosols have a short lifetime so if future generations stop injecting them, this is what will happen. The animation shows that the surface temperature increased by global warming is recovered after the removal of aerosols.

To run this simulation, we first initialise the gas with the largest optical depth parameters that will be encountered in the simulation, indicated by the final suffix here. This is so the pressure grid has enough points in areas where there is a high mass concentration. We then set the optical depth parameters to their start values and update the optical depth grid values using the function update_grid.

We next find the analytical equilibrium solution for the starting optical depth and use this as the starting point for the simulation. In the data dictionary, we also include the optical depth because it is changing, so we want to add it to the animation.

When running the simulation, we call take_time_step(t, T_eqb.copy(), changing_tau) during times where the optical depth is changing. This is because we must include the fact that changing_tau=True so the function recomputes the optical depths at each time. At times when we just want to reach equilibrium, we just run evolve_to_equilibrium, providing the current state of the data dictionary. We stop the simulation either after 10 years have passed or we have reached the final equilibrium after the aerosols have been removed.

Example 2 - increasing $\tau_{lw}$ for thermosphere

This second example shows the result of increasing the surface $\tau_{lw}$ from 4 to 8 for the thermosphere example. Basically, the simulation achieves equilibrium with the current $\tau$ parameters and then increases tau_lw_surface by 0.1 and repeats.

The starting profile is the profile that the previous thermosphere simulation ended on and this is the temperature profile shown in orange on the left plot. I think the initial increase in mesopause temperature is because the convergence threshold of delta_net_flux_thresh = 1e-3 is much stricter here. After this start though, you can see a trend for surface temperature to increase, mesopause to narrow and its temperature to decrease and tropopause to shift up.

The narrowing occurs because there is very little interaction of long radiation with gas above $\tau_{lw}=1$ so to match the boundary conditions at the top of the atmosphere, the radiation coming from the $\tau_{lw}=1$ level must be constant and thus so must the temperature according to $\sigma_B T^4$. So when the $\tau_{lw}=1$ level moves to a lowe pressure, the temperature profile has to shift up aswell. This consequently pushes the tropopause closer to the mesopause.

In other words, as optical depth increases, the highest pressure value where the blue line intersects the red line moves to a lower pressure. To intersect at this reduced pressure level, the surface temperature has to increase. I wonder whether it could somehow jump though so that the emitting pressure level is the intersection in the stratosphere instead of in the troposphere, in which case it could be achieved without a change in surface temperature.

I tried to investigate this by, increasing the long wave gas concentration in the upper atmosphere in the cell after this. There are some big fluctuations in temperature so I am not sure if it is reliable but it shows that there is now only one intersection between the red and blue lines and this is very high in the atmosphere. So now it is not the temperature profile near the surface that is important to satisfy the top of atmosphere boundary condition but the profile near the top of the atmosphere.

Latitudinal dependence

You can also choose to include a latitudinal dependence of the incoming solar radiation by setting ny>1. This will then make ny latitudes evenly spaced between $-90^{\circ}$ and $+90^{\circ}$. The flux that arrives at each latitude is determined from the annually averaged value as given by equation 13.2.6 in Atmospheric Circulation dynamics and General Circulation Models. The latitude columns do not interact and the default initial condition is so that each latitude column is isothermal and the incoming short wave radiation is equal to the outgoing long wave radiation as the top of the atmosphere. You can also specify the albedo to be different for different latitudes. At the moment, it assumes that the optical depth distribution is the same for each latitude column.

The animation shows the temperature as a function of pressure and latitude as well as the surface temperature as a function of latitude. You can also provide the optical depth but the flux cannot be plotted as there are too many values when considering each latitude.

The below example just repeats the aerosols in the stratosphere simulation but with 30 latitude bins.